/* Definition of quadrature formulas
   Jan Stebel
*/

#ifndef _QUADRATURE_H__
#define _QUADRATURE_H__



/**
  Structure for storing the quadrature formulas
*/
template<unsigned int dim>
struct Quadrature
{
  int n;
  const double *p[dim], *w;
};

// 3D Gaussian quadrature formulas for tetrahedra (exact for 4-th order polynomials)
const double Q3d11px[] = { 0.2500000, 0.7857143, 0.07142857, 0.07142857, 0.07142857, 0.1005964, 0.1005964, 0.1005964, 0.3994034, 0.3994034, 0.3994034, 0, 1, 0, 0, 0.5, 0.5, 0,   0,   0,   0.5 };
const double Q3d11py[] = { 0.2500000, 0.07142857, 0.7857143, 0.07142857, 0.07142857, 0.1005964, 0.3994034, 0.3994034, 0.1005964, 0.1005964, 0.3994034, 0, 0, 1, 0, 0,   0.5, 0.5, 0,   0.5, 0   };
const double Q3d11pz[] = { 0.2500000, 0.07142857, 0.07142857, 0.7857143, 0.07142857, 0.3994034, 0.1005964, 0.3994034, 0.1005964, 0.3994034, 0.1005964, 0, 0, 0, 1, 0,   0,   0,   0.5, 0.5, 0.5 };
const double Q3d11pw[] = { -0.01315556, 0.007622222, 0.007622222, 0.007622222, 0.007622222, 0.02488889, 0.02488889, 0.02488889, 0.02488889, 0.02488889, 0.02488889, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
const Quadrature<3> Q3d11p = { 11, { Q3d11px, Q3d11py, Q3d11pz }, Q3d11pw };

const double Q3d15px[] = { 0.2500000, 0.0000000, 0.3333333, 0.3333333, 0.3333333, 0.7272727, 0.09090909, 0.09090909, 0.09090909, 0.4334498, 0.4334498, 0.4334498, 0.06655015, 0.06655015, 0.06655015, 0, 1, 0, 0, 0.5, 0.5, 0,   0,   0,   0.5  };
const double Q3d15py[] = { 0.2500000, 0.3333333, 0.0000000, 0.3333333, 0.3333333, 0.09090909, 0.7272727, 0.09090909, 0.09090909, 0.4334498, 0.06655015, 0.06655015, 0.4334498, 0.4334498, 0.06655015, 0, 0, 1, 0, 0,   0.5, 0.5, 0,   0.5, 0    };
const double Q3d15pz[] = { 0.2500000, 0.3333333, 0.3333333, 0.0000000, 0.3333333, 0.09090909, 0.09090909, 0.7272727, 0.09090909, 0.06655015, 0.4334498, 0.06655015, 0.4334498, 0.06655015, 0.4334498, 0, 0, 0, 1, 0,   0,   0,   0.5, 0.5, 0.5  };
const double Q3d15pw[] = { 0.03028368, 0.006026786, 0.006026786, 0.006026786, 0.006026786, 0.01164525, 0.01164525, 0.01164525, 0.01164525, 0.01094914, 0.01094914, 0.01094914, 0.01094914, 0.01094914, 0.01094914, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
const Quadrature<3> Q3d15p = { 15, { Q3d15px, Q3d15py, Q3d15pz }, Q3d15pw };

const double Q3d45px[] = { 0.2500000,  0.6175872,  0.1274709,  0.1274709,  0.1274709,
                           0.9037635,  0.03207883, 0.03207883, 0.03207883, 0.4502229,
			   0.4502229,  0.4502229,  0.04977710, 0.04977710, 0.04977710,
			   0.3162696,  0.3162696,  0.3162696,  0.1837304,  0.1837304,
			   0.1837304,  0.5132800,  0.5132800,  0.5132800,  0.02291779,
			   0.02291779, 0.02291779, 0.2319011,  0.2319011,  0.2319011,
			   0.2319011,  0.2319011,  0.2319011,  0.1937465,  0.1937465,
			   0.1937465,  0.7303134,  0.7303134,  0.7303134,  0.03797005,
			   0.03797005, 0.03797005, 0.03797005, 0.03797005, 0.03797005,
			   0, 1, 0, 0, 0.5, 0.5, 0,   0,   0,   0.5 };
const double Q3d45py[] = { 0.2500000,  0.1274709,  0.6175872,  0.1274709,  0.1274709,
                           0.03207883, 0.9037635,  0.03207883, 0.03207883, 0.4502229,
			   0.04977710, 0.04977710, 0.4502229,  0.4502229,  0.04977710,
			   0.3162696,  0.1837304,  0.1837304,  0.3162696,  0.3162696,
			   0.1837304,  0.02291779, 0.2319011,  0.2319011,  0.5132800,
			   0.2319011,  0.2319011,  0.5132800,  0.5132800,  0.02291779,
			   0.02291779, 0.2319011,  0.2319011,  0.7303134,  0.03797005,
			   0.03797005, 0.1937465,  0.03797005, 0.03797005, 0.1937465,
			   0.1937465,  0.7303134,  0.7303134,  0.03797005, 0.03797005,
			   0, 0, 1, 0, 0,   0.5, 0.5, 0,   0.5, 0 };
const double Q3d45pz[] = { 0.2500000,  0.1274709,  0.1274709,  0.6175872,  0.1274709,
                           0.03207883, 0.03207883, 0.9037635,  0.03207883, 0.04977710,
			   0.4502229,  0.04977710, 0.4502229,  0.04977710, 0.4502229,
			   0.1837304,  0.3162696,  0.1837304,  0.3162696,  0.1837304,
			   0.3162696,  0.2319011,  0.02291779, 0.2319011,  0.2319011,
			   0.5132800,  0.2319011,  0.02291779, 0.2319011,  0.5132800, 
			   0.2319011,  0.5132800,  0.02291779, 0.03797005, 0.7303134,
			   0.03797005, 0.03797005, 0.1937465,  0.03797005, 0.7303134,
			   0.03797005, 0.1937465,  0.03797005, 0.1937465,  0.7303134,
			   0, 0, 0, 1, 0,   0,   0,   0.5, 0.5, 0.5 };
const double Q3d45pw[] = { -0.03932701,   0.004081316,  0.004081316,  0.004081316,  0.004081316,
                            0.0006580868, 0.0006580868, 0.0006580868, 0.0006580868, 0.004384259,
                            0.004384259,  0.004384259,  0.004384259,  0.004384259,  0.004384259,
                            0.01383006,   0.01383006,   0.01383006,   0.01383006,   0.01383006,
                            0.01383006,   0.004240437,  0.004240437,  0.004240437,  0.004240437,
			    0.004240437,  0.004240437,  0.004240437,  0.004240437,  0.004240437,
			    0.004240437,  0.004240437,  0.004240437,  0.002238740,  0.002238740,
			    0.002238740,  0.002238740,  0.002238740,  0.002238740,  0.002238740,
			    0.002238740,  0.002238740,  0.002238740,  0.002238740,  0.002238740,
			    1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
const Quadrature<3> Q3d45p = { 45, { Q3d45px, Q3d45py, Q3d45pz }, Q3d45pw };

// 2D Gaussian quadrature formulas for triangles
const double Q2d4px[] = { 1./3,    0.2,   0.2,   0.6,   0, 1, 0, 0.5, 0.5, 0 };
const double Q2d4py[] = { 1./3,    0.6,   0.2,   0.2,   0, 0, 1, 0,   0.5, 0.5 };
const double Q2d4pw[] = { -27./96, 25./96, 25./96, 25./96, 1, 1, 1, 1,   1,   1 };
const Quadrature<2> Q2d4p = { 4, { Q2d4px, Q2d4py }, Q2d4pw };

const double Q2d6px[] = { 0.8168476,  0.09157621, 0.09157621, 0.1081030, 0.4459485, 0.4459485, 0, 1, 0, 0.5, 0.5, 0 };
const double Q2d6py[] = { 0.09157621, 0.8168476 , 0.09157621, 0.4459485, 0.1081030, 0.4459485, 0, 0, 1, 0,   0.5, 0.5 };
const double Q2d6pw[] = { 0.05497587, 0.05497587, 0.05497587, 0.1116908, 0.1116908, 0.1116908, 1, 1, 1, 1, 1, 1 };
const Quadrature<2> Q2d6p = { 6, { Q2d6px, Q2d6py }, Q2d6pw };

const double Q2d7px[] = { 1.012865e-01, 1.012865e-01, 7.974270e-01, 4.701421e-01, 4.701421e-01, 5.971587e-02, 3.333333e-01, 0, 1, 0, 0.5, 0.5, 0 };
const double Q2d7py[] = { 1.012865e-01, 7.974270e-01, 1.012865e-01, 4.701421e-01, 5.971587e-02, 4.701421e-01, 3.333333e-01, 0, 0, 1, 0,   0.5, 0.5 };
const double Q2d7pw[] = { 6.296959e-02, 6.296959e-02, 6.296959e-02, 6.619708e-02, 6.619708e-02, 6.619708e-02, 1.125000e-01, 1, 1, 1, 1,   1,   1 };
const Quadrature<2> Q2d7p = { 7, { Q2d7px, Q2d7py }, Q2d7pw };

// 1D Gaussian quadrature formulas
const double Q1d2px[] = { -sqrt(3)/3, sqrt(3)/3, -1., 1., 0. };
const double Q1d2pw[] = { 1., 1., 1., 1., 1. };
const Quadrature<1> Q1d2p = { 2, { Q1d2px }, Q1d2pw };

const double Q1d3px[] = { 0., -sqrt(15)/5, sqrt(15)/5, -1., 1., 0. };
const double Q1d3pw[] = { 8./9, 5./9, 5./9, 1., 1., 1. };
const Quadrature<1> Q1d3p = { 3, { Q1d3px }, Q1d3pw };

const double Q1d4px[] = { -sqrt(525.-70.*sqrt(30.))/35., sqrt(525.-70.*sqrt(30.))/35., -sqrt(525.+70.*sqrt(30.))/35., sqrt(525.+70.*sqrt(30.))/35., -1., 1., 0. };
const double Q1d4pw[] = { (18.+sqrt(30.))/36., (18.+sqrt(30.))/36., (18.-sqrt(30.))/36., (18.-sqrt(30.))/36., 1., 1., 1. };
const Quadrature<1> Q1d4p = { 4, { Q1d4px }, Q1d4pw };

const double Q1d10px[] = { -0.1488743389816312108848260, 0.1488743389816312108848260,
                           -0.4333953941292471907992659, 0.4333953941292471907992659,
                           -0.6794095682990244062343274, 0.6794095682990244062343274,
                           -0.8650633666889845107320967, 0.8650633666889845107320967,
                           -0.9739065285171717200779640, 0.9739065285171717200779640,
                           -1., 1., 0. };
const double Q1d10pw[] = { 0.2955242247147528701738930, 0.2955242247147528701738930,
                           0.2692667193099963550912269, 0.2692667193099963550912269,
                           0.2190863625159820439955349, 0.2190863625159820439955349,
                           0.1494513491505805931457763, 0.1494513491505805931457763,
                           0.0666713443086881375935688, 0.0666713443086881375935688,
                           1, 1, 1 };
const Quadrature<1> Q1d10p = { 10, { Q1d10px }, Q1d10pw };










#endif
